#ifndef _EBFV1_ELLIPTIC_EQUATION_H_
#define _EBFV1_ELLIPTIC_EQUATION_H_

#include "Elliptic_equation.h"
//#include "TMesh.h"

namespace PRS{
	// User defined Matrix-vector product for matrix-free operations
	int MatMultUser(Mat mat, Vec u, Vec y);

	/**
	 * For EBFV1 elliptic formulation, a matrix-free procedure is applied to solve pressure field. A pointer-function is used to tell Petsc how to make
	 * a matrix-vector product and all matrices and vector are passed to this function as a structure.
	 */
	struct MatVec_struct{
		Mat G;				// global matrix (all mesh nodes)
		Mat *E;				// global matrix (all mesh nodes)
		Mat *F;				// global matrix (all mesh nodes)
		Mat EF;				// global matrix (all mesh nodes)
		Mat G_freeNodes;	// global matrix (free mesh nodes only)
		Mat *E_freeNodes;	// global matrix (free mesh nodes only)
		Mat *F_freeNodes;	// global matrix (free mesh nodes only)
		Mat EF_freeNodes;	// global matrix (free mesh nodes only)
		Vec RHS;
		Vec solution;		// solution vector
		int ndom;
		int nrows;
		int *rows;

		Vec z;
		int F_nrows;
		int F_ncols;
	};

	class EBFV1_elliptic : public Elliptic_equation{
	public:
		EBFV1_elliptic();
		EBFV1_elliptic(PhysicPropData *,SimulatorParameters *,MeshData *);
		~EBFV1_elliptic();
		double solver(TMesh* pTMesh_data);

	private:
		bool allocateMatVec;					// check if matrices have been assembled
		Mat matrix;								// PETSc matrix for matrix-free procedure
		MatVec_struct *pMatVec;					// Pointer for matrices and vector structure
		int PETSc_assemblyMatrix(Mat);			// After fill the matrix, it must be assembled (PETSC)
		int assembly_EFG_freeNodes(TMesh *pTMesh_data);
		double solver_MatFree();
		int zeroMatricesAndVectors();

		/*!brief Creates an array of indices (row = freeNodesIndices[i]) for direct access to pTMesh_data->vertices[dom].setValue(row,0,data).
		 *       This minimizes data time transfer from PETSc vector 'pMatVec->solution' to VertexData pointers.
		 *
		 */
		int createAuxiliaryMappingForSolution(TMesh *pTMesh_data);

		// Fill matrices E,G and F for a specific edge and domain. Theses function are called inside a loop of domains
		void calculateAuxiliaryDataForMatrices_EFG(TMesh *pTMesh_data);

		/*
		 * Matrix G is equivalent to the matrix of discretization of elliptic terms of DFVC method in a dual Voronoi mesh. G is assembled with a loop over all
		 * edges. In parallel, if an edge is on the partition boundary, Gij must be divided by (NRC) Number of Remote Copies. Each partition will add Gij/NRC
		 * to the global matrix E and the resulting Gij contribution will be the same as it was assembled in serial.
		 */
		int divergence_G(Mat G, int id0, int id1, double MobIJ, double aux);

		/*
		 *  Matrix E is responsible to project the gradient over an orthogonal plane to edge. E is assembled with a loop over all edges. In parallel, if an edge
		 *  is on the partition boundary, Eij must be divided by (NRC) Number of Remote Copies. Each partition will add Eij/NRC to the global matrix E and the
		 *  resulting Eij contribution will be the same as it was assembled in serial.
		 */
		int divergence_E(Mat G, int dim, int id0, int id1, double MobIJ, EdgeData* edata);
		int gradient_F_omega(Mat F, int dim, int id0, int id1, EdgeData *edata);
		int gradient_F_gamma(Mat F, TMesh* pTMesh_data);
		double assemblyRHS(TMesh* pTMesh_data);
		int setRHS_Neumann(TMesh* pTMesh_data, Vec RHS);
		int setRHS_Dirichlet(Mat A, Vec RHS);
		double assemblyMatrices(TMesh *pTMesh_data);
		int allocateMatricesAndVectors(TMesh* pTMesh_data);
		double updatePressure(TMesh* pTMesh_data);
		double pressureGradient(TMesh* pTMesh_data);

		// MatMultUser is called several time by Petsc during KSP solver until convergence be reached. Multiplication => A*u = y, where A = EF+G
		static int MatMultUser(Mat mat, Vec u, Vec y){
			void *ctx;
			PetscErrorCode ierr = MatShellGetContext(mat,&ctx);  CHKERRQ(ierr);
			MatVec_struct* mats = (MatVec_struct*)(ctx);
			ierr = MatMult(mats->G_freeNodes,u,y); CHKERRQ(ierr);
			for (int i=0; i<mats->ndom; i++){
				ierr = MatMultAdd(mats->EF_freeNodes,u,y,y); CHKERRQ(ierr);
			}
			return ierr;
		}


		PhysicPropData *pPPData;
		SimulatorParameters *pSimPar;
		MeshData *pMData;
		Matrix<int> freeNodesIndices;
	};
}
#endif
